#################################################################################
# Legacies of Collective Action: Why Families of Disappeared Mobilized in Mexico#
# Laura Lopez-Perez                                                             #
# University of Notre Dame                                                      #
# Violence and Transitional Justice Lab                                         #
# Eliminating Violence Against Women Lab                                        #
# llopez7@nd.edu                                                                #
# Last update: 08/05/24                                                         #
#################################################################################

rm(list = ls())
colectivos <- read.csv("colectivos.csv")

###### Plot - timeseries ####

library(ggplot2)
library(dplyr)
library(date)
library(lubridate)

colectivosyear$year <- as.numeric(as.character(colectivosyear$year))
colectivosyear$colectivos <- as.numeric(as.character(colectivosyear$colectivos))

fecha <- colectivosyear %>%
  mutate(fecha = make_datetime(year))

colectivosyear$date <- as.Date(fecha$fecha)

p <- ggplot(fecha, aes(x=year, y=colectivos_a)) +
  geom_line() +
  geom_point() +
  xlab("") +
  ylab("Colectivos of families of disappeared persons (total)") +
  geom_vline(xintercept=2006, linetype="longdash") +
  geom_vline(xintercept=2011, linetype="longdash") +
  geom_vline(xintercept=2014, linetype="longdash") +
  theme_minimal()
p + theme(legend.title = element_blank()) +
  scale_x_continuous(breaks=seq(2001, 2020, 1)) +
  scale_y_continuous(breaks=seq(0, 160, 10))


###### Plot - timeseries #####

fecha <- colectivos %>% 
  select(year) %>% 
  mutate(fecha = make_datetime(year))

colectivos$date <- as.Date(fecha$fecha)

df1 <- colectivos %>%
  filter(state %in% c("Coahuila", "Tamaulipas", "Zacatecas", "Baja California", "Guerrero"))

p <- ggplot(df1, aes(x=year, y=accumulated, group=state)) +
  geom_line(aes(color=state)) +
  xlab("") +
  ylab("Disappeared persons (accumulated rate)") +
  geom_vline(xintercept=2001, linetype="longdash") +
  geom_vline(xintercept=2007, linetype="longdash") +
  geom_vline(xintercept=2009, linetype="longdash") +
  geom_vline(xintercept=2011, linetype="longdash") +
  geom_vline(xintercept=2014, linetype="longdash") +
  theme_minimal()
p + theme(legend.title = element_blank()) +
  scale_x_continuous(breaks=seq(2001, 2014, 1)) +
  scale_y_continuous(breaks=seq(0, 5500, 500))

#### Create variable of national and local elections####
colectivos$elections <- 0
colectivos$elections[colectivos$year==2006 | colectivos$year==2012] <- 1

colectivos$stateel <- 0
colectivos$stateel[colectivos$year==2001 & colectivos$state=="Baja California" | 
                     colectivos$state=="Michoac?n" | colectivos$state=="Yucat?n"] <- 1
colectivos$stateel[colectivos$year==2003 & colectivos$state=="Campeche" |
                     colectivos$state=="Colima" | colectivos$state=="Guanajuato" |
                     colectivos$state=="Nuevo Leon" | colectivos$state=="Quer?taro" |
                     colectivos$state=="San Luis Potos?" | colectivos$state=="Sonora"] <- 1
colectivos$stateel[colectivos$year==2004 & colectivos$state=="Aguascalientes" | 
                     colectivos$state=="Chihuahua" | colectivos$state=="Durango" |
                     colectivos$state=="Puebla" | colectivos$state=="Oaxaca" | 
                     colectivos$state=="Sinaloa" | colectivos$state=="Tamaulipas" |
                     colectivos$state=="Tlaxcala" | colectivos$state=="Veracruz" |
                     colectivos$state=="Zacatecas"] <- 1
colectivos$stateel[colectivos$year==2005 & colectivos$state=="Baja California Sur" |
                     colectivos$state=="Coahuila" | colectivos$state=="Colima" |
                     colectivos$state=="Guerrero" | colectivos$state=="Hidalgo" |
                     colectivos$state=="Estado de M?xico" | colectivos$state=="Nayarit" |
                     colectivos$state=="Quintana Roo"] <- 1
colectivos$stateel[colectivos$year==2006 & colectivos$state=="Chiapas" |
                     colectivos$state=="Ciudad de M?xico" | colectivos$state=="Guanajuato" |
                     colectivos$state=="Jalisco" | colectivos$state=="Morelos" |
                     colectivos$state=="Tabasco"] <- 1
colectivos$stateel[colectivos$year==2007 & colectivos$state=="Baja California" |
                     colectivos$state=="Michoac?n" | colectivos$state=="Yucat?n"] <- 1
colectivos$stateel[colectivos$year==2009 & colectivos$state=="Campeche" |
                     colectivos$state=="Colima" | colectivos$state=="Nuevo Le?n" |
                     colectivos$state=="Quer?taro" | colectivos$state=="San Luis Potos?" |
                     colectivos$state=="Sonora"] <- 1
colectivos$stateel[colectivos$year==2010 & colectivos$state=="Aguascalientes" | 
                     colectivos$state=="Chihuahua" | colectivos$state=="Durango" |
                     colectivos$state=="Puebla" | colectivos$state=="Oaxaca" |
                     colectivos$state=="Sinaloa" | colectivos$state=="Tamaulipas" |
                     colectivos$state=="Tlaxcala" | colectivos$state=="Veracruz" |
                     colectivos$state=="Zacatecas"] <- 1
colectivos$stateel[colectivos$year==2011 & colectivos$state=="Baja California Sur" |
                     colectivos$state=="Coahuila" | colectivos$state=="Guerrero" |
                     colectivos$state=="Hidalgo" | colectivos$state=="Estado de M?xico" |
                     colectivos$state=="Nayarit" | colectivos$state=="Quintana Roo"] <- 1
colectivos$stateel[colectivos$year==2012 & colectivos$state=="Chiapas" |
                     colectivos$state=="Ciudad de M?xico" | colectivos$state=="Guanajuato" |
                     colectivos$state=="Jalisco" | colectivos$state=="Morelos" |
                     colectivos$state=="Tabasco"] <- 1
colectivos$stateel[colectivos$year==2013 & colectivos$state=="Michoac?n"] <- 1
colectivos$stateel[colectivos$year==2014 & colectivos$state=="Baja California"] <- 1

################### Logit models #########
m1 <- glm(colectivos$colectivos ~ colectivos$hr_org, family=binomial('logit'))
m2 <- glm(colectivos$colectivos ~ colectivos$hr_org + 
            colectivos$disaprate + colectivos$acumulated_rate + colectivos$homiciderate,
          family=binomial('logit'))
m3 <- glm(colectivos$colectivos ~ colectivos$hr_org + 
            colectivos$disaprate + colectivos$acumulated_rate + colectivos$homiciderate +
            colectivos$income_osc + colectivos$elections + colectivos$stateel,
            family=binomial('logit'))

library(stargazer)
stargazer(m1, m2, m3, type='text')


########## Cross-sectional model on the effect of the MPJD #######

## Subset dataset for cross-sectional analysis
mpjd <- colectivos[colectivos$year>=2011 & colectivos$year<=2013,]

csm1 <- glm(mpjd$colectivos ~ mpjd$mpjd,
            family=binomial('logit'))
csm2 <- glm(mpjd$colectivos ~ mpjd$mpjd + mpjd$hr_org +
              mpjd$disaprate + mpjd$acumulated_rate,
            family=binomial('logit'))
csm3 <- glm(mpjd$colectivos ~ mpjd$mpjd + mpjd$hr_org +
              mpjd$disaprate + mpjd$acumulated_rate + 
              mpjd$income_osc + mpjd$elections + mpjd$stateel,
            family=binomial('logit'))

library(stargazer)
stargazer(csm1, csm2, csm3, type='text')

############################ Instrumental variable analysis #########################
# Generate subsections of the dataset
iv1 <- colectivos[colectivos$year<=2010, ]
iv2 <- colectivos[colectivos$year>=2011, ]

fitX.LZ <- glm(formula=hr_org~liberationtheology + femospp, family="poisson", data=colectivos)
fitY.LX <- glm(formula=colectivos~hr_org + femospp, family="binomial", data=colectivos)
fitIV.1 <- ivglm(estmethod="ts", fitX.LZ=fitX.LZ, fitY.LX=fitY.LX, data=colectivos, 
               ctrl=TRUE) 
summary(fitIV.1)

fitX.LZ1 <- glm(formula=hr_org~liberationtheology + femospp, family="poisson", data=iv1)
fitY.LX1 <- glm(formula=colectivos~hr_org + femospp, family="binomial", data=iv1)
fitIV1.1 <- ivglm(estmethod="ts", fitX.LZ=fitX.LZ1, fitY.LX=fitY.LX1, data=iv1, 
                ctrl=TRUE) 
summary(fitIV1.1)

fitX.LZ2 <- glm(formula=hr_org~liberationtheology + femospp, family="poisson", data=iv2)
fitY.LX2 <- glm(formula=colectivos~hr_org + femospp, family="binomial", data=iv2)
fitIV2.1 <- ivglm(estmethod="ts", fitX.LZ=fitX.LZ2, fitY.LX=fitY.LX2, data=iv2, 
                ctrl=TRUE) 
summary(fitIV2.1)


########## Instrumental Variable - Assumptions Tests #######

## Instrument strength
t1 <-  t.test(colectivos$comites, colectivos$hr_org)
cov(colectivos$comites, colectivos$hr_org)
t2 <- t.test(colectivos$liberationtheology, colectivos$hr_org)
cov(colectivos$liberationtheology, colectivos$hr_org)

## Graph confidence intervals
#Instrument strenght
library("ggpubr")
ggscatter(iv1, x = "comites", y = "hr_org", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Comités of families of forcefully disappeared persons", 
          ylab = "Local human  rights organizations")

ggscatter(iv2, x = "liberationtheology", y = "hr_org", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Liberation Theology Movement", 
          ylab = "Local human  rights organizations")


## Exclusion restriction
# Comités
t1 <- t.test(colectivos$comites, colectivos$colectivos, conf.level = 0.95)

mean_CI <-(t1$conf.int[2]+t1$conf.int[1])/2
d <- data.frame(x = 0, 
                y = (t1$conf.int[2] + t1$conf.int[1]) / 2,
                ymin = t1$conf.int[1], 
                ymax = t1$conf.int[2])

plot1 <- ggplot(d)+ 
  geom_errorbar(aes(x, ymin = ymin, ymax = ymax), width = .1) +
  geom_point(aes(x, y), size = 4)+
  geom_hline(yintercept = 0, linetype = 3, size =1.5) +
  xlim(-.2, .2) +
  ylim(-0.8,0.8) +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())


# Liberation Theology
t3 <- t.test(iv1$liberationtheology, colectivos$comites, iv1$colectivos, conf.level = 0.95)
t4 <- t.test(iv2$liberationtheology, colectivos$comites, iv2$colectivos, conf.level = 0.95)

#Graph confidence intervals
mean_CI <-(t3$conf.int[2]+t3$conf.int[1])/2
d <- data.frame(x = 0, 
                y = (t3$conf.int[2] + t3$conf.int[1]) / 2,
                ymin = t3$conf.int[1], 
                ymax = t3$conf.int[2])

plot3 <- ggplot(d)+ 
  geom_errorbar(aes(x, ymin = ymin, ymax = ymax), width = .1) +
  geom_point(aes(x, y), size = 4)+
  geom_hline(yintercept = 0, linetype = 3, size =1.5) +
  xlim(-.2, .2) +
  ylim(-0.8,0.8) +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

mean_CI<-(t4$conf.int[2]+t4$conf.int[1])/2
d <- data.frame(x = 0, 
                y = (t4$conf.int[2] + t4$conf.int[1]) / 2,
                ymin = t4$conf.int[1], 
                ymax = t4$conf.int[2])

plot4 <- ggplot(d)+ 
  geom_errorbar(aes(x, ymin = ymin, ymax = ymax), width = .1) +
  geom_point(aes(x, y), size = 4)+
  geom_hline(yintercept = 0, linetype = 3, size =1.5) +
  xlim(-.2, .2) +
  ylim(-0.8,0.8) +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())


